A mean-field kinetic lattice gas model of electrochemical cells 
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We develop Electrochemical Mean-Field Kinetic Equations (EMFKE) to simulate electrochemical 
cells. We start from a microscopic lattice-gas model with charged particles, and build mean-field 
kinetic equations following the lines of earlier work for neutral particles. We include the Poisson 
equation to account for the influence of the electric field on ion migration, and oxido-reduction 
processes on the electrode surfaces to allow for growth and dissolution. We confirm the viability 
of our approach by simulating (i) the electrochemical equilibrium at flat electrodes, which displays 
the correct charged double-layer, (ii) the growth kinetics of one-dimensional electrochemical cells 
during growth and dissolution, and (iii) electrochemical dendrites in two dimensions. 

PACS numbers: 82.45.Qr, 05.70.Ln, 64.60.-i 



I. INTRODUCTION 

Electrochemical phenomena are ubiquitous in nature 
and technology. They play a fundamental role in many 
materials science problems of high practical relevance, 
such as corrosion, electrodeposition of parts and circuitry, 
and battery technology pj . Electron transfer is also cru- 
cially involved in many biochemical reactions 0]. For 
physicists, electrodeposition is interesting because it can 
lead to the spontaneous formation of highly complex self- 
organized patterns. Such branched aggregates that are 
often found on the surfaces of natural minerals have fas- 
cinated scientists for centuries because of their plant- like 
appearance [3j- 

Under well-controlled laboratory conditions, a vari- 
ety of different patterns can be generated, ranging from 
compact crystals to highly branched dendritic aggregates 
that can be fractals (see Figure or densely branched. 
While dendrites rarely appear in traditional applications 
of electrodeposition such as metalization and electroplat- 
ing that take place close to equilibrium, they can be of 
practical importance far from equilibrium, for example in 
battery technology 4] : Figure|21shows the growth of den- 
drites in a lithium battery. Since dendrites can perforate 
insulating layers and lead to shortcuts when they reach 
the counterelectrode, they must be eliminated as much as 
possible. To achieve this goal, a good understanding of 
the various phenomena occurring in electrochemical cells 
is necessary. The purpose of our paper is to develop a 
microscopic model for electrodeposition that can be used 
to elucidate some of the aspects of dendrite formation. 

To fix the ideas, let us consider the simplest possible 
situation that corresponds to the experimental setup of 
Figs. PI] and [2] two electrodes made of the same metal 
are plunged in an electrolyte containing ions of the same 
metal. When a potential difference is applied by an ex- 
ternal generator, as a response an ionic current flows 
through the electrolyte. Close to the cathode, ions are 
reduced by electron transfer from the electrode and form 
a growing deposit. The inverse takes place at the an- 
ode: metal is dissolved, and new ions are formed. For 



FIG. 1: Electrodeposition of copper on a glass substrate 
(reprinted with permission from V. Fleury). The quasi- 
two-dimensional cell is made of two copper electrodes and 
a CuS04 solution (sample size 3x2 mm 2 ). Only a small part 
of the sample is shown. 



liquid electrolytes (such as used in Fig. G), the inhomo- 
geneities in the ion concentrations lead to strong convec- 
tive flows [HIES HI- In contrast, for gel-like electrolytes 
(such as used in Fig.|2J), convection is suppressed and 
the important ingredients for the description of the ion 
transport in the electrolyte are diffusion and migration 
of the ions. At the metal-electrolyte interface, there ex- 
ists in general a charged double layer that is much thicker 
than the microscopic solid-electrolyte interface, but much 
smaller than a typical cell dimension. There is a consid- 
erable amount of work on electrochemical cells based on 
a macroscopic viewpoint [lfj in which the microscopic 
interface is replaced by a mathematically sharp surface. 
Such models are generally unable to handle the geometri- 
cal complexity of a fully developed dendrite. In addition, 
this approach has to rely on phenomenological models to 
incorporate reaction kinetics at the interfaces. 

In the statistical physics community, considerable in- 
terest in electrodeposition was spurned by the fact 
that the geometrical structure of fractal electrodeposits 
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FIG. 2: Dendritic growth at the interface between a lithium 
electrode and a polymer electrolyte (Claire Brissot, PhD the- 
sis, Ecole Polytechnique, Palaiseau, France, 1998) 



ated by the diffusion-limited aggregation (DLA) model 
|15| . in which random walkers irreversibly stick to the 
growing aggregate. Subsequently, many refinements of 
the DLA model were developed to generate various pat- 
terns, for example by the inclusion of surface tension [16j, 
anisotropic growth rules and noise reduction |l7l ITsI ] , in- 
troduction of a uniform drift to mimic an electric field 
[Til I20L I2H ]. and combinations of discrete aggregation 
and continuous convection models to study the influence 
of fluid convection j^. However, these models greatly 
simplify the underlying physics (for example, the charged 
double layers are not included) and can hence not yield 
detailed information on the relation between growth con- 
ditions and characteristics of the growth structures such 
as growth speed, branch thickness, and overall structure. 
The same is true for a recent mean-field model that does 
not explicitly contain the charged ions |23j . 

Our goal is to build a simplified microscopic model 
for electrodeposition that includes enough of the salient 
physics to make contact with the macroscopic view of 
electrochemistry. The electrochemical mean-field ki- 
netic equations (EMFKE) developed here are based on a 
lattice-gas model with simple microscopic evolution rules 
that contains charged particles, coupled to a discretized 
version of the Poisson equation. Lattice gas models have 
been used previously in the context of electrochemistry 
to simulate phenomena located on the electrode surfaces, 
such as adsorption or underpotential deposition [24| , and 
for studies of ionic transport at liquid-liquid interfaces 
|25|: however, there exists, to our knowledge, no theo- 
retical study of the behavior of an entire electrochemical 
cell based on a microscopic model. 

To investigate the dynamics of the lattice model, we 
extend the formalism of Mean-Field-Kinetic-Equations 
(MFKE) IHIillHl that has been used to study numer- 
ous transport and growth phenomena in alloys, including 
diffusion and ordering kinetics 29], spinodal decomposi- 
tion H(], EI an d dendritic growth |32|: some preliminary 
results on the extension to electrochemistry have been 
published in fiefs. (3^, The key feature of this ap- 
proach is that the microscopic particle currents can be 
written in the mean-field approximation as the product 



of a mobility times the gradient of chemical potentials, 
the latter being the appropriate thermodynamic driving 
forces. The formalism can be generalized in a natural 
way to charged particles. The driving forces for particle 
currents are then the gradients of the electrochemical po- 
tential. As a consequence, the resulting model displays 
the correct electrochemical equilibrium at the interfaces. 
While it obviously still contains strong simplifications, 
it is much closer to the basic microscopic physics and 
chemistry than the DLA-type models and allows us to 
establish a direct link between a microscopic model and 
the well-established macroscopic phenomenological equa- 
tions. Quite remarkably, the equations of motion of our 
model also share many common features with a phase- 
field formulation of the electrochemical interface that was 
developed very recently 35]. This indicates that mean- 
field equations of simplified models can also be useful to 
understand the connection between phase-field models 
and a more microscopic viewpoint. 

In this paper, we will first present the lattice gas 
model and derive the electrochemical mean-field ki- 
netic equations (Section 2). In section 3, we carry 
out one-dimensional simulations to demonstrate that the 
model leads to the correct equilibrium at the electrode- 
electrolyte interface. We also calculate one-dimensional 
steady state solutions for moving interfaces that are in 
good agreement with the macroscopic continuous model 
of Chazalviel |36j| . Finally, we show some preliminary 
simulations of dendritic growth in two-dimensional cells. 
Section 4 contains a brief discussion and conclusion. 



II. MODEL 
A. Lattice gas model 

We consider an electrochemical cell made of a di- 
lute binary electrolyte and two metallic electrodes of the 
same metal. No supporting electrolyte is included in the 
present study. These conditions correspond to the ex- 
perimental situation in Figs. ^ and [21 (more precisely, to 
Fig. [3 since we do not include convection in our model) . 
The electrodes are modeled by a lattice that reflects the 
underlying crystalline structure, and whose sites are oc- 
cupied by metallic atoms or vacancies. It is convenient to 
represent the electrolyte by the same lattice, occupied by 
a solvent, cations, anions, or vacancies. Although there 
is no physical lattice in the liquid, its presence here does 
not play a role due to the high dilution of the ions. 

For simplicity, we consider here a two-dimensional lat- 
tice gas on a square lattice with lattice spacing a (see 
Figure 13 ■ The cations M + give metallic atoms M° after 
reduction, while the anions A~ are supposed to be non 
electroactive. The solvent S is neutral, but can interact 
through short-range interactions with the other species 
and with itself. We specify a microscopic configuration 
by the set {n} of the occupation numbers n£: on each site 
k: ng = 1 if k is occupied by species a = M°, M + ,A~ , S 
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or a vacancy v, and otherwise. We suppose steric ex- 
clusion between the different species, that is, a given site 
can be occupied by only one species or it can be empty 
(vacancy): 

a 

Finally, as will become clear below, it is convenient to 
introduce electrons in the metallic electrodes. They have 
a particular status that will be discussed later. 

Two very different types of interactions have to be 
considered. The short-range interactions (including, 
for example, van der Waals forces, solvatation effects, 
and chemical interactions) are modeled here by nearest- 
neighbor interactions e a @ between species a and j3 (with 
the convention that a positive e a/3 corresponds to an at- 
tractive interaction) . Interaction energies with vacancies 
are taken to be zero. To take into account the long- 
range electrostatic interaction, it would be possible to 
introduce appropriate interaction energies with farther 
neighbors; however, this procedure becomes cumbersome 
with increasing interaction range. Instead, we consider 
a "coarse-grained" electrical potential Vk defined on the 
lattice sites k (that is, a potential that has been smoothed 
over distances smaller than a) and that obeys Poisson's 
equation with the simplest nearest-neighbor discretiza- 
tion for the Laplace operator, 

a a— + 

where e is the dielectric constant (for simplicity, we take 
a constant that is the same for all species), q a is the 
electric charge of species a, and d the spatial dimension. 
This equation has to be solved in the electrolyte, that is, 
outside of the metal clusters connected to the ends of the 
cell, and subject to the boundary condition of constant 
potential in each electrode. 

The configurations evolve when one particle (metal, 
ion, or solvent) jumps to one of its vacant nearest neigh- 
bor sites. In principle, we could also include an exchange 
process between occupied nearest neighbor sites. This 
leads to more complicated kinetic equations and has not 
been considered here. To specify the jump rates, we 
assume that the atoms perform activated jumps. The 
height of the activation barrier depends on the local bind- 
ing energy, that is, the number and type of bonds that 
need to be broken (see Fig. 0J, and, for charged particles, 
on the local electric field that shifts the barrier height as 
shown in Fig. [S] The result for the jump rate from site 
k to site k + a is 

S£ )k+ a(W) = ^P^EE^+a^ 

X eX P (£f( V * - Vk+aOV ( 3 ) 




FIG. 3: The lattice gas model. A fixed potential difference is 
applied across the cell. The ions in the electrolyte are sub- 
mitted to an electric field Ek (and hence a force Fk = f/Ek) 
at their lattice site position k. The various species have short 
range interactions (here, attractive interactions are consid- 
ered between solvent and ions, solvent and solvent, and metal 
and metal). Electron transfer takes place on the electrode 
surfaces. 




FIG. 4: In the lattice gas model, an A atom makes activated 
jumps to empty nearest neighbor sites. The barrier it has to 
overcome depends on its interactions (white links) with its 
nearest neighbor atoms. 



U k -U k+a /2 q(V k -V k+a )/2 qE k+a/2 .a/2 




X 



k k+a 



FIG. 5: In presence of an electric field E, a potential energy 
which varies, to first order, like gE.x along the jump path 
x, is superimposed to the local potential seen by a moving 
particle (with charge q). 
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where kT is the (fixed) thermal energy, and w a is a fixed 
jump frequency that may be different for each species. 
Here and in the following, symbols with tildes will de- 
note electrochemical quantities that are sensitive to the 
electric potential. 



a+a' 



P a' 



+ 5 |1(VW-VW) 



(8) 



B. Mean-field kinetic equations 

The derivation of the Electrochemical Mean-Field Ki- 
netic Equations (EMFKE) follows the same procedure 
as for neutral particles [2|| ■ We first write the 

Boolean kinetic equations on the lattice, starting from 
the general master equation, 

§-P({n},t) = £[W({n'}-{n})P({n'},t) 
{»'} 

-W{{n}^{n'})P{{n},t)] (4) 

which gives the probability P to find a given configu- 
ration {n} at time t. W({n} — > {n 1 }) is the rate of 
evolution from configuration {n} to configuration {n'}. 
We are interested in the time evolution of the average 
concentrations of the various species (a = +, — , 0, S, v), 



£n£P({n},t). 

{n} 



(5) 



In the absence of electrochemical processes at the elec- 
trodes, the number of each type of particles remains con- 
stant, and the kinetic equation of the average concentra- 
tion has the structure of a conservation equation, 



2_J ^k,k+a' 



with the current of species a on link {k, k - 
by 



(6) 



a} defined 



k,k4 



k.k4 



>(W> 



k"k+a 



"k+a.k 



({n})» 



k+a 



(7) 

with w£ k+a ({n}) given by Eq. ©. The factor n^nj^ +a 
(and n£ +a n£ for the reverse jump) means that for a jump 
of a to be possible, the start site must be occupied by 
species a, while the target site must be empty (occupied 
by a vacancy v). 

In the mean-field approximation, the occupation num- 
bers n£ in the above expressions are replaced by their av- 
erage p£ . This replacement is not unique because of dif- 
ferent possible choices for the factorization of the occupa- 
tion number operators j2?J • A convenient choice |23, H2 
is the direct replacement of all occupation numbers by 
their averages in Eqs. (7J and © , which leads to 



J k,k+a 



PkPk + aexp(- ^EE e " 



^k+a' 



2kT 



13 a' 



(Vk - V k+a 



This expression can be rewritten without any further ap- 
proximation as the product of an electrochemical bond 
mobility My times the (discrete) gradient of an electro- 
chemical potential Jlf, 



Ji 



k,k+a 



i®aMk: 



(9) 



where S a is a difference operator acting on the site coor- 
dinates, 2) a /k = /k+a — /k- The electrochemical poten- 
tial, 



Mk = Mk+<Ai 



<7^ k ,(10) 



is the sum of three contributions, a local energy due to 
the interaction of species a with its local environment, an 
entropy term (these two constitute the chemical potential 
/ig), and an electrostatic energy q a V\ i . The presence of 
the vacancy concentration in the denominator of the en- 
tropic contributions comes from the constraint of Eq. Q . 
The mobility along a bond i -j is given by 



JIJ k,k+a 



^PkPk+aexp — 



-she 



SaMk 

2kT 



(11) 



where we have used the notation she u = smhu/u (close 
to equilibrium, pt£+a = and she [J) a ^/2fcT] = 1 ). 

In a dilute electrolyte, where the concentration of 
species a — +, — , 0, and v is low, neglecting all the terms 
of order larger than two in the concentrations and dis- 
crete gradients, the current J^ k+a simplifies, 

(12) 

where < = w§ exp (-(1/fcT) E/3 E a ' ^Vk+a') • This 
is the discrete form of the continuous macroscopic ex- 
pression 



-D aa grad c Q - D^grad c v 



D aa —cE, 
kT ' 



(13) 



with a concentration-dependent diffusion coefficient 
D aa = a 2 w a , an off-diagonal diffusion coefficient D av = 
a 2 w a c a /c v associated with the gradient of the vacancy 
concentration, and the correspondence between a d- 
dimensional cubic lattice and the d— dimensional continu- 



ous space Jk,k+s 



J,Pk 



c. When the vacancy 



concentration is homogeneous, the off-diagonal term is 
absent, and Eq. (|I3|) is exactly the classical continuous 
description of ion diffusion and migration. 
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So far, the mean-field equations are quite similar to 
the ones governing the evolution of multi-component al- 
loys. However, two important new elements have now to 
be taken into account: (i) how to calculate the electric 
potential in the mean-field model, and (ii) the electro- 
chemical reactions at the electrodes. 



C. The Poisson problem 

The electric potential has to be obtained, as before, by 
resolving Poisson's equation. However, some additional 
considerations are necessary, because it turns out that 
the treatment of the boundary conditions at the elec- 
trode surfaces becomes non-trivial in the mean-field con- 
text. To understand this, it is useful to start with some 
comments on the consequences of the mean-field approx- 
imation. The parameters of the model that control the 
phase diagram are the various interaction energies and 
the temperature. Obviously, we want to create and main- 
tain two distinct phases, namely the metallic electrodes 
and the electrolyte. Hence, the interaction energies need 
to be chosen such that phase separation occurs. As usual 
for this kind of lattice models, there exists a critical tem- 
perature for phase separation. Close below the critical 
point, the equilibrium concentrations of the two phases 
are close to each other, that is, there are many metal 
particles in the electrolyte phase, and vice versa. The 
concentration of these minority species decreases with 
temperature. In the original lattice gas model, it is then 
possible for low enough temperatures to identify the ge- 
ometry of the bulk phases with the connected clusters of 
metal and electrolyte species. 

In the mean-field approximation, the concentration 
variables (occupation probabilities) are continuous, and 
each species has a non-vanishing concentration at each 
site. This means that the electrode always contains small 
quantities of solvent and ions, and the electrolyte con- 
tains a small quantity of metal. Of course, these concen- 
trations can be made arbitrarily small by lowering the 
temperature. However, we then face another difficulty. 
The interface between the two bulk phases, which was 
the sharp boundary between connected clusters in the 
discrete model, is now diffuse with a characteristic thick- 
ness that depends on the temperature. For low temper- 
atures, this thickness becomes smaller than the lattice 
spacing, which leads to strong lattice effects on static 
and dynamic properties of the interface [3?1 l38l l3^| : the 
surface tension and mobility of the interface depend on 
its position and orientation with respect to the lattice, 
and for very low temperatures the interface can be en- 
tirely pinned in certain directions. As a consequence, we 
have to make a compromise in choosing the temperature: 
it must be low enough to obtain reasonably low concen- 
trations of the minority species, but high enough to avoid 
lattice pinning and to obtain a reasonably diffuse inter- 
face (that extends over a few lattice sites). 

This has consequences for the resolution of Poisson's 



equation. The boundary conditions of constant electric 
potential in the electrodes are easy to impose in a model 
with sharp boundaries. For diffuse interfaces, we have to 
decide where the electrodes end. Another way to state 
the problem is to remark that the existence of an elec- 
tric field in the electrolyte creates surface charges in the 
conductor that are localized at the surface. In a system 
with diffuse interfaces, this surface charge is "smeared 
out" over the thickness of the interface, and we need a 
method to determine this charge distribution. 

We solve these problems by the introduction of elec- 
trons that are free to move within the metallic electrodes 
and solve the Poisson equation for all the charges, in- 
cluding electrons. More precisely, we denote by p£ the 
deviation from the neutral state expressed in electrons 
per site. Hence, p£ > corresponds to an excess of elec- 
trons, p£ < to an electron deficit. Their chemical po- 
tential is defined by 

K = E F + q ^ + ^- y (14) 

T>{Ep) being the density of electronic states at the Fermi 
level Ep of the metal, and q e = —e. This corresponds to 
the screening approximation of Thomas- Fermi [40] . The 
electronic current is then written as an electronic mobil- 
ity times the discrete gradient of the chemical potential, 

^k,k+a = -A?k,k+aSa^ k , (15) 

and the time evolution of the excess electron concentra- 
tion is 

a 

Note that the choice of this dynamics is somewhat ar- 
bitrary; it could be replaced by any other process that 
leads to an equilibrium of the electron distribution with 
the ionic charge distribution on a time scale that is much 
faster than the characteristic time of evolution of the 
other species. 

One important feature is that the electrons must re- 
main in the electrodes; otherwise, a "leakage" current 
may occur. We have tested two different possibilities to 
achieve this goal. We can either consider that the elec- 
trons are in a potential well in the metallic regions such 
that they remain confined in the metal, or we can suppose 
that their mobility vanishes outside the metallic regions. 
Here, we have chosen the second possibility and written 
the mobility in the form, 

— in^ 

^k,k+a = j^f(j>l)f(pl+J, (17) 

where w e is a constant frequency prefactor, and / is an 
interpolation function that is equal to 1 for large metal 
concentrations and falls to zero for low metal concentra- 
tions. With this choice, the electronic jump probability 
is important only if nearest neighbor sites k and k + a 
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have a large enough probability to be occupied by metal- 
lic atoms. We have used for / 



fip) 



tanh[(p - ffc)/£] + tanh[p c /£] 
tanh[(l - p c )/£] + tanh[p c /£] ' 



(18) 



a monotonous function that varies from when p = to 
1 for p = 1, with a rapid increase through an interval in 
p of order £ centered around some concentration p c that 
is reminiscent of a percolation threshold. This interpo- 
lation is motivated by the fact that the metallic region 
must be dense enough to be connected and thus allow 
the electrons to propagate. 

To determine the electrostatic potential, we solve the 
mean-field version of Poisson's equation, including the 
new contribution from the electrons, 



(19) 



a— +. — ,e 



This equation is now solved in the whole system, includ- 
ing the electrodes. Note that we still use a constant per- 
mittivity e; however, the phase-dependent mobility of the 
electrons makes the resistivity in the electrolyte much 
higher than in the metal. 

This method provides a fast way to calculate the sur- 
face charges on the electrodes. As will be shown below, it 
works perfectly well at equilibrium. However, in out-of- 
equilibrium simulations, a problem appears on the side 
of the anode where the metal is dissolved. Since the 
mobility rapidly decreases with the metal concentration, 
electrons present on the metallic site before dissolution 
may be trapped in the electrolyte, leading to spurious 
electronic charges in the bulk. We have solved this prob- 
lem in a phenomenological way by adding a term to the 
evolution equation for the electrons that relaxes the elec- 
tronic charge to zero in the electrolyte, 



dpi 
dt 



= -E J k,k + a-^ e [l-/,(PkM, (20) 



where f r (j>°) is the same interpolation function as /, but 
with different parameters £ r and p c ,r- With a convenient 
choice of these parameters, "electron relaxation" occurs 
only in the liquid. 



species. Reduction of cations on a site k of the cathode 
appears when 



Mk + Mk+a > Mk> 



(22) 



otherwise, the metal is oxidized. Consequently, we define 
Ck,k+a as the current of electronic charges from k + a to 
k (current of positive charges from k to k + a ) reducing 
the cations on site k (resp. electronic current issued from 
the oxidation of the metal) via a corresponding elimina- 
tion (resp. creation) of electrons on site k + a. Following 
llOll. we can write the reaction rate, 



Ck,k+a - w k,k-] 



exp 



K 



kT 



exp # 



(23) 



This corresponds to an activated electronic charge trans- 
fer between the metal surface and the nearest neighboring 
cation. The total reduction rate on site k is the sum of 
all the reaction paths J2 a "k.k+a- The coefficient u>£ k+a 
can be determined by comparison with the mesoscopic 
theory of Butler- Volmcr (see the Appendix) . 

The currents Ck.k+a have non negligible values only 
on the interfaces. However, as discussed before, there 
are small concentrations of ions in the metal and metal 
atoms in the electrolyte; this would lead to undesirable 
contributions of the reaction term inside the bulk phases. 
Therefore, we use the same interpolation function f(p) as 
for the electron mobility and write 



w 



k,k+a 



™*(! -/(*£)) /(Pk+a), 



(24) 



where w* is again a constant frequency factor. In this 
way, the transfer is localized around the metal-electrolyte 
interface, where occupied metallic sites and electrolyte 
sites are neighbors. To completely suppress the electro- 
chemical processes in the bulk phases, we set the transfer 
current to zero when the product of the two interpolation 
prefactors is smaller than 10~ 4 . 



E. Summary 

After all the different pieces are combined, the com- 
plete set of equations is 



dp+ 



D. Electron transfer 



Q t ~ J k,k+a / / 



Ok, k+a, 



(25) 



On the electrode surfaces, electron transfer takes 
places. More precisely, metallic cations M + located in 
the electrolyte may receive an electron from a neighbor- 
ing metallic site and be reduced; in turn, metal atoms in 
contact with the electrolyte may reject an electron to a 
neighboring metallic site and become an ion, 



M+(k) + e~(k + a) ^ M°(k) 



(21) 



The direction of the transfer depends on the relative mag- 
nitude of the electrochemical potentials of the involved 



M 

dt 



xi j k.k+ a + x 



Ck,k+a, 



dPk 
dt 

dpi 

dt 



X/ ^k,k+a' 



k,k+a' 



(26) 



(27) 



(28) 
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-jjf = ~Y1 J k.k+a - W e [l - f r (pi)]Pk ~ °k+a,k, 
a a 

(29) 

5>k +a -4T/ k = -^ «■ (30) 

a a— + ,— ,e 

The local concentration of the species M° and M + is 
modified by transport (diffusion and migration in the 
electric field) and by the chemical reaction; for the other 
species A~ and S, only transport is present. The reac- 
tion term is active only in the solid-electrolyte interface 
and has of course an opposite sign in Eqs. Ij25(l and (|26J) . 
and no contributions in Eqs. Ij27(l and (|28|l . 

Equations (|25|) through l|29[l are integrated in time by a 
simple Euler scheme (with constant or variable timestep). 
Two methods are used to solve Poisson's equation. Equa- 
tion (H3) can either be solved for each timestep using a 
conjugate gradient method, or it can be converted into 
a diffusion equation with a source term and solved by 
simple relaxation; the second method generally is com- 
putationally faster and converges well for interfaces that 
move slowly enough. 

III. NUMERICAL TESTS 
A. Chemical equilibrium 

The complete EMFKE model involves a considerable 
number of parameters. In particular, we need to make a 
choice for all the short-range interaction energies. They 
determine the phase diagram of the four-component sys- 
tem that consists of metal, solvent, anions, and cations. 
To obtain meaningful results, the model has to be used 
not too far from a two-phase equilibrium between a 
metal-rich electrode and a liquid phase rich in solvent 
and ions. For an arbitrary choice of interaction ener- 
gies, obtaining the equilibrium concentrations for all of 
the species in the two phases is not a trivial task. The 
conditions for a two-phase equilibrium are (i) equal elec- 
trochemical potentials for each independent species (four 
in our case), and (ii) equal grand potential. This yields 
five conditions for eight unknowns (four equilibrium con- 
centrations in each phase); consequently, three degrees 
of freedom remain that may be fixed, for example the 
concentrations of ions and metal in the liquid. This con- 
stitutes a set of five coupled nonlinear equations for the 
remaining unknowns. 

To simplify the task of finding an equilibrium configu- 
ration, we start from a simpler system, namely a mixture 
of metal, solvent, and vacancies, that is equivalent to a bi- 
nary alloy with vacancies (ABv). For a symmetric inter- 
action matrix, e 00 = e ss , the phase diagram can be deter- 
mined analytically |3l| . A two-phase equilibrium exists 
with p s Uq = p° sol and pf iq — p s sol . The equilibrium concen- 
trations can be expressed in terms of the new variables 

P = Ptiq+Pliq = P S sol +P Q sol and Q = Pliq-Pliq = P2oJ -Psol 



and the reduced interaction energy e = (e°°+e ss — 2e 0s )/2 
as 

P = Q/tanh(jgV (31) 

where z is the coordination number [z = 4 for a square 
lattice in two dimensions). For the rest of the section, 
we choose e 00 /kT = e ss /kT = 1, e 0s = 0, and Q = 0.9, 
which gives pi iq = p° sol = 0.92528, p% q = p* sol = 0.02528, 
and pl q =p» sol = .0.04944. 

Next, this equilibrium is perturbed by the addition of 
ions. If the interaction energies are taken identical for 
both species of ions, no separation of charge occurs at 
equilibrium and in the absence of an applied voltage, and 
we can for the moment omit the Poisson equation. We 
take e +s /kT = e~ s /kT = 1 and e++ = e~ = e+° = 
= 0. This means that the ions are attracted by the 
solvent, and that the energy density of the liquid does not 
change upon addition of ions. As long as the equilibria 
of the main components are not appreciably modified, 
we expect an equilibrium distribution of the ions that 
satisfies 

We tested this prediction by performing simulations of 
one-dimensional half cells of length 40, in which 20 lattice 
sites of solid are in contact with 20 lattice sites of liquid. 
In each phase, all concentration values are initially set 
to the expected equilibrium values. On the solid side, 
we apply no-flux boundary conditions, whereas at the 
liquid side the concentrations are kept constant. This 
corresponds to a finite-size isolated electrode plunged in 
an infinite solution. The electrochemical reaction is sup- 
pressed by setting w* — 0; all other frequency factors w a 
are set to 1. The equations are integrated with a fixed 
timestep of At = 1. 

Two ways of adding ions were tested. In the first case, 
p\ iq is kept constant and a small percentage of solvent 
is replaced by ions. Within a few 10 6 timesteps, the 
chemical potential differences across the cell fall below 
10 -9 and the evolution virtually stops. The concentra- 
tion profiles obtained for pf iq = pj iq = 0.01 are shown 
in Fig. |SJ In Fig. \7\ we plot the ion concentration in 
the solid versus the ion concentration in the liquid. For 
low concentrations, Eq. Ij32fl is well satisfied. For concen- 
trations larger than « 10 -3 , deviations appear because 
the concentrations of metal, solvent and vacancies are 
shifted. However, even for pf iq = 0.01 this deviation is 
less than 5%. 

In the second case, ions were added in replacement 
of vacancies {p s Hq and p® iq were kept constant). Here, 
the solvent-metal equilibrium is shifted more rapidly due 
to the strong dependence of the chemical potentials on 
the vacancy concentration. For ion concentrations larger 
than 0.003, the interface never relaxes, but starts to grow 
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FIG. 6: Concentration profiles obtained from the relaxation 
of a step profile (20 sites of solid in contact with 20 sites of 
electrolyte) with ion concentration in the liquid oip^ iq = 0.01. 
The other parameters are given in the text. 




ion concentration (liquid) 



FIG. 7: Ion concentration in the solid p sol in equilibrium with 
a liquid of fixed ion content p^ iq . 



by incorporating metal that is transported to the inter- 
face from the solution by chemical diffusion. Of course, 
the fact that the system does not reach equilibrium is due 
to the fact that we have fixed one degree of freedom too 
much, namely all four concentrations in the liquid. Con- 
sequently, even the solutions found by the first method 
are not, strictly speaking, equilibrium solutions; however, 
they are sufficiently close to equilibrium for all practical 
purposes and will hence be used as initial conditions in 
the following. 



FIG. 8: Concentration profiles with the same parameters as 
in Fig. |S| but with a potential difference AV = — 0.5fcT/e 
between metal and solution. 



B. Poisson equation and screening: perfectly 
polarizable electrode 

To test our method for creating interface charges, we 
started from an equilibrated half cell as calculated in 
the preceding subsection, added a potential difference be- 
tween electrode and solution, and opened the metal side 
of the cell for a current of electrons (but not of other 
species). Since the electron transfer frequency was kept 
to zero (w* = 0), and hence no electrochemical processes 
can take place, this situation corresponds to a grounded 
perfectly polarizable electrode. 

For simplicity, we have restricted our attention here 
to the case of monovalent ions (q + — e, q~ = — e). In 
our tests, we want the thickness of the charged double 
layer to be a few lattice constants. This can be achieved 
by choosing the dielectric constant e = 0.05e 2 / (akT). 
Furthermore, we choose for the density of states at the 
Fermi level V(E F ) = 1000/fcT. Since kTV{E F ) > 1, 
the electrochemical potential of the electrons depends 
only weakly on the surface charge. The parameters of 
the interpolation functions are fixed to p c = 0.5, £ = 0.1 
for the electron mobility (this leads to an electron mo- 
bility that is 8 orders of magnitude smaller in the liq- 
uid than in the electrodes) and p c . r = 0.03, £ r = 0.001 
for the electron relaxation (see Eq. (|20|l : this corre- 
sponds to an interpolation function with a sharp increase 
at a concentration slightly larger than the equilibrium 
metal concentration in the liquid). The resulting inter- 
face profile for pf iq — 0.01 and a potential difference of 
AV = V meta i - Vu q md = -0.5kT/e is plotted in Fig. 
As expected, the cations are attracted to the electrode, 
whereas the anions are repelled. 

In Fig. [51 we show the charge distribution that cor- 
responds to the profile of Fig. |SJ Several features are 
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FIG. 9: Charge distribution (q = q a 
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noteworthy. The small homogeneous ion concentrations 
inside the electrodes are now different for positive and 
negative ions. This is a direct consequence of the global 
equilibrium: since the electrochemical potentials are con- 
stant throughout the whole system and the ion concen- 
trations in the liquid are fixed, the ion concentrations 
in the solid are multiplied by factors exp(— q a AV/kT) 
with respect to the reference state AV = 0. The result- 
ing "background charge" in the bulk electrode is exactly 
compensated by the electrons that have flown in from the 
metal side of the system. Since the resulting total charge 
inside the bulk electrode is zero, the potential is constant, 
which was the original purpose of introducing the elec- 
trons. Inside the diffuse interface between electrode and 
solution, we find a surface charge that is located almost 
exclusively on one site. 

In the liquid, we find the diffuse charged layer that 
is predicted by the classic sharp-interface model of Gouy 
and Chapman. For an electrolyte that has singly charged 
ions, q + = —q~ = e, the potential and the deviations of 
the concentrations from the bulk liquid values all decay 
exponentially with the distance from the interface. For 
example, V(x) = AV exp[— n(x — Xi n t)] with a decay con- 
stant 



ekT ' 



(33) 



the inverse of the Debye screening length, and Xi nt the 
position of the (sharp) interface. In Fig. we show 
the potential profile across a cell subjected to a potential 
difference of AV — O.lfcT/e for different ion concentra- 
tions. We obtain values for k by an exponential fit of 
V(x), using Xint (the position of the "sharp interface" 
extrapolated from the far field) as an adjustable param- 
eter. The results are in good agreement with the predic- 
tions of Eq. The potential follows an exponential 
law over several orders of magnitude. However, when 
the potential approaches its reference value in the liquid, 



FIG. 10: Potential V as a function of x for AV = O.lfeT/e 
and pf iq =0.0001 (solid line), pf iq =0.001 (dashed line), and 
p^ 9 =0.01 (dash-dotted line). The thick line segments are 
exponential fits (V(x) = AVexp[— n(x — !;„()]) that yield 
k = 0.0626, 0.199, and 0.622 for the inverse of the Debye 
length. The predictions of Eq. £g) are 0.0632, 0.200, and 
0.632, respectively. The sizes of the simulated cells were 200, 
100, and 50 lattice spacings. 



deviations from the exponential behavior occur. This is 
due to the fact that the buildup of the ion boundary 
layer at the interface also modifies the chemical poten- 
tials for solvent and metal by small amounts. Since the 
electrochemical potential of the ions depends on the sol- 
vent and vacancy concentrations, a complete equilibrium 
solution must contain these effects. Close to the inter- 
face, they are negligible in comparison to the contribu- 
tion of the electric potential; only far from the interface, 
when V — Vu qu id becomes of the same order as the small 
shifts in the chemical potentials, the deviations from an 
exponential become appreciable. Since this occurs for 
fractions of the potential drop that are less than 10~ 3 of 
AV, these effects can safely be neglected in the further 
analysis. 

The resulting total charge distribution has the ex- 
pected double-layer structure, with a much sharper decay 
in the electrode than in the liquid. We have checked that 
the integral of the total charge through the interface is 
zero to numerical precision. 



C. Electrochemical equilibrium: Nernst law 

The electrochemical equilibrium between electrodes 
and electrolyte depends on the choice of the Fermi en- 
ergy Ep . Indeed, consider a chemical equilibrium in the 
absence of an electric potential such as the one shown in 
Fig. El In this state, the electrochemical potential of the 
ions is constant across the cell. If we fix Ep to be exactly 
equal to the difference of the chemical potentials of metal 
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and cations, 

E% = - fi+> (34) 

the reaction currents a are strictly zero everywhere as 
long as the potential remains constant. Therefore, the 
original profile remains an equilibrium state for arbitrary 
transfer frequency w* . In contrast, if Ep is set to a value 
different from E F , the electron transfer starts and drives 
the system to a new equilibrium that exhibits a potential 
difference. 

The final potential difference will be close to Ep — E F . 
This can be deduced as follows. For fixed concentrations 
and electric potential in the liquid, the electrochemical 
potentials of the ions and the metal are fixed. Therefore, 
the only way to achieve equilibrium is to adjust the elec- 
tric potential in the electrode such that the two terms 
in Eq. I|23l) balance. Now consider the electrochemical 
potential of the electrons. The excess electronic charge 
p e is always small, \p e \ <C eAV/kT, even at the sur- 
face, as can be seen in Fig. [5] Since, in addition, we 
have chosen T){Ep)kT ^> 1, Jl e is only weakly dependent 
on the local electron density. Neglecting the last term 
in Eq. 1141 the electron transfer rate becomes zero when 
E F - eAV = E° F . 

We have checked this prediction by performing simu- 
lations in isolated cells: the concentrations of all species 
and the potential were kept fixed at the liquid side, 
and the system was closed at the solid side for parti- 
cles and electrons, thus enforcing zero electrical current. 
For pf lq = 0.01, Eq. (JSJl yields E%/kT = A.U75U. For 
Ep = E F , the system remained at constant potential, 
whereas for Ep ^ E F , the predicted potential difference 
developed up to an error of less than 10~ 5 . Since all of 
the charges can be removed from the system by applying 
an external potential difference that is just the negative 
of (Ep—E F ) I e, the choice of Ep determines the potential 
of zero charge (PZC) of the electrochemical interface. 

Of course, the value of E F as defined above depends 
on the ion concentration in the liquid. Therefore, for 
fixed Ep, there is a well-defined ion concentration p^ e f 
for which the potential difference is zero. When the con- 
centration in the liquid is varied, the potential difference 
follows Nernst's law, 

AV=— ln^-. (35) 

e Pref 

In Fig. EU we plot AV vs pf lq for E° F /kT = 4.447514 
(corresponding to pp. — 0.01). The results of our model 
are in perfect agreement with Eq. (|35|l . 

D. Growth 

To investigate growth and dissolution, we use one- 
dimensional cells with a metallic electrode at both ends. 
The system is closed on both sides for all species (ions, 




ion concentration (liquid) 



FIG. 11: Equilibrium potential difference through an iso- 
lated cell versus ion concentration in the liquid for fixed 
E F = 4.447514 kT. The solid line corresponds to Eq. 

metal, solvent) but open for electrons. To drive the in- 
terfaces, we apply a fixed potential difference AV across 
the cell. For eAV/kT ranging from 1 to 10 and a cell of 
size 100 with two electrodes of thickness 10, the system 
reaches a steady state within 5 x 10 6 timesteps (At = 1). 
Within 2 x 10 7 timesteps, the interface advances between 
3 and 15 sites, depending on the voltage. We measure the 
ionic current in the center of the cell by time-averaging 
over the second half of the runs. 

Given the rapid variations of ion concentrations and 
electronic charges through the interface (see for example 
Fig. IHL one could expect strong lattice pinning effects 
[371 l38t Isfjj . We have checked that quantities such as 
the electronic surface charge, the total transfer rate (that 
is, the transfer currents <7k,k+a summed up through the 
interface) and the ion current indeed do vary as the in- 
terface advances through the lattice. However, for the 
parameters chosen here, the amplitudes of these lattice 
oscillations never exceeded a few percent. 

The choice of the various time constants does not in- 
fluence the final results for equilibrium states. In con- 
trast, for growth simulations, they have to be fixed in 
order to achieve the desired physical conditions. In par- 
ticular, this is true for the various jump rates and the 
electron transfer frequency. The electric current in the 
electrolyte is carried by the mobile ions. As mentioned 
above, the buildup of the ionic boundary layers shifts 
the chemical potential of metal, solvent, and vacancies 
as well. Since, in our model, the concentration of metal 
in the electrolyte is comparable to the ion concentrations, 
this leads to neutral diffusion currents that are unphys- 
ical. To lower their magnitude, the jump frequency for 
the metal was chosen much smaller than for the ions: 
w + = w~ = 1, w° = 10~ 3 . 

In Fig. 1121 we plot the ionic current in the center of 
the cell as a function of the transfer frequency w* at 
fixed voltage of 10 kT/e. For low transfer frequencies, the 
growth is limited by the reaction kinetics at the interface, 
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FIG. 12: Ion current in the center of the cell as a function 
of the electron transfer frequency w* for fixed voltage AV = 
WkT/e. 
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FIG. 13: Ion current in the center of the cell as a function of 
the driving voltage for fixed w* = 10 -4 . 



FIG. 14: Concentration profiles (top) and potential profile 
(bottom) across a 100-site cell with initial ion concentration 
of 0.01 subjected to a potential difference of A V = kT/e. The 
snapshot was taken at t = 2 x 10 7 . 




and the current strongly depends on the value of w* . For 
increasing w* , the growth becomes limited by transport 
in the bulk, and the current is almost independent of w* . 
We are mostly interested in the latter regime. Therefore, 
we fix in the following w* = 10~ 4 . 

The current-voltage curve for our model cell is shown 
in Fig. It is strongly nonlinear. Indeed, two very 
different regimes are covered by these simulations. This 
can be appreciated when looking at the ion and potential 
profiles in Figs, d and US For AV = kT/e (Fig. 113$!. the 
ion concentrations remain of the same order of magnitude 
as the initial concentration (p = 0.01), and except for 
the two double layers close to the interfaces, the liquid is 
neutral. An important part of the potential drop occurs 
in these double layers; in the bulk electrolyte, the poten- 
tial profile is smooth and almost linear. In contrast, for 
AV = 10fcT/e (Fig.UHJ), the neighborhood of the cathode 



FIG. 15: Same as Fig. 1141 but for a potential difference of 
AV = 10 kT/e. 



has been completely depleted of the anions (p~ ~ 10~ 5 
at the cathode), and a charged zone extends well beyond 
the thickness of the equilibrium double layers. Most of 
the potential drop occurs close to the cathode. Since the 
conductivity in the space-charge zone is lowered due to 
the low ion concentration, a considerable increase in the 
overall voltage AV leads only to a moderate increase of 
the current, as seen in Fig. EH All of these findings are 
in good agreement with the macroscopic one-dimensional 
calculations of Chazalviel [3y|. According to this work, 
the extended space charge is crucial for the emergence 
of ramified growth: the strong electric field close to the 
surface leads to an instability of the flat front, and one- 
dimensional calculations become invalid. 
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E. Two-dimensional simulations: dendritic growth 



We present now an example for a preliminary simu- 
lation of a two-dimensional sample. The purpose is to 
show that our model can indeed lead to the emergence 
of dendritic structures; however, the conditions that we 
can simulate are far from typical experimental situations. 
The main reason is that in experiments, the branches of 
ramified aggregates have typically a thickness in the mi- 
cron range, whereas the lattice of our model represents 
a crystal lattice with spacing of the order Angstroms. 
It is clear that huge simulation cells would be necessary 
to observe instabilities and ramified growth on realistic 
scales. To obtain a computationally tractable problem, 
we have to work with unrealistically high driving forces, 
since this is known to reduce the characteristic scales of 
branched growth structures. Therefore, we use a much 
higher driving potential than in the previous simulations, 
namely 100 kT/e. This corresponds to about 2.5 V at 
room temperature, which is a fairly typical value; how- 
ever, this potential difference is applied through a cell 
that is, as before, 100 lattice sites long, which corre- 
sponds to a length of a few nanometers. Therefore, the 
electric fields are much higher in our simulation than in 
reality. 

Under these conditions, the behavior of the moving in- 
terface is quite sensitive to the model parameters, and in 
particular to the frequency factors for the electron trans- 
fer, w*, and the metal jump frequency w° (as before, we 
take w + = ur = 1 as a reference value). To assure that 
the growth is still transport-limited for the higher driv- 
ing force, w* has to be chosen large enough; for reaction- 
limited growth, no morphological instability occurs. The 
metal jump frequency has to be chosen carefully. On 
the one hand, if it is too high, bumps on the surface 
are smoothed out too rapidly by surface diffusion and/or 
an evaporation-condensation mechanism. On the other 
hand, if it is too low, the diffusion of metal on the solid 
side of the interface becomes so slow that the interface 
profile cannot be maintained, and the metal grows at a 
concentration far below its equilibrium value [41|. 

In Fig. 1161 we show an example computed with w* = 
6 x 10~ 3 , w° = 10~ 3 , and w e = 10~ 3 . The cell has a size 
40 x 100 lattice sites, with periodic boundary conditions 
parallel to the interfaces. The simulation was started 
from a flat interface, with random shifts of the metal 
concentration in the interfaces to trigger the instability. 
When one layer of the anode was dissolved, the whole 
cell was shifted backward by one site in order to keep the 
electrolyte in the center of the cell. It can be seen that a 
bump grows on the interface and develops into a finger- 
like structure. Other bumps that initially develop on the 
interface are screened. The whole region that surrounds 
the dendrite is depleted of ions. An extended charged 
region forms ahead of the tip. When the tip gets closer to 
the anode, the electric field increases and leads to growth 
of the metal at unphysically low concentrations. 




FIG. 16: Snapshots of the evolution of a 100 x 40 electrochem- 
ical cell; the parameters are given in the text. The electrodes 
are in black; the cathode is at the bottom and grows up- 
ward (the cell is recentered during the simulation). In the 
electrolyte, the grayscale indicates the concentration of ions, 
with white corresponding to the maximal ion concentration. 
The white "contour" at the electrode surfaces is the charged 
double layer; the gray region in front of the dendrite has been 
almost completely depleted of ions. 



IV. CONCLUSION 

In summary, we have shown here that, start- 
ing from a simple microscopic model, it is possible 
to build Electrochemical-Mean-Field-Kinetic Equations 
(EMFKE) that are able to reproduce qualitatively the 
behavior of electrochemical cells. Both the charged dou- 
ble layers present at equilibrium and the extended space 
charge that develops during growth are correctly repro- 
duced. Dendritic structures can be simulated, albeit for 
unrealistic parameters. Hence, the EMFKE contain the 
fundamental ingredients that are necessary to simulate 
dendritic growth by electrodeposition. 

Our model shares many common features with a recent 
phase-field formulation of electrodeposition 35]. The 
phase-field method, originally developed in the context 
of solidification in the 1980s p3. l43ll4^ |. is a continuous 
model of phase transitions that uses an auxiliary indica- 
tor field, the phase field, to distinguish between the dif- 
ferent thermodynamic phases (here, electrodes and elec- 
trolyte). A phenomenological equation of motion for the 
phase field is usually derived from a free energy func- 
tional. In our EMFKE approach, the role of the phase 
field is played by the metal concentration. 

Since the phase-field method is phenomenological, it 
has a certain freedom of choice for the dynamics of the 
phase field itself, and can hence avoid the problems that 
arise in our model due to the presence of metal in the elec- 
trolyte. However, no direct link to a microscopic model is 
established, which is the strength of our approach. Our 
equations still contain some phenomenological elements, 
in particular, the interpolation function for electron dif- 
fusivity and reaction rates. A more realistic modeling of 
the processes involving electrons is needed to overcome 
this limitation. With this perspective, our approach may 
constitute a useful link between microscopic models and 
phase-field models. 
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Reaction coordinate 



for the redox reaction, and (if we suppose that the fre- 
quency prefactors are the same) the corresponding rates 
are 



k rec i — A exp — 



AG re( z(</>) 



A exp 



AG ox f» 



For some potential 
AG red (0 o ) = AG 



kT 1 kT 

(A2) 

, the two barriers are equal, 
and hence k rer i = k nT such that 



the total reaction current j is zero. When <j) ^ ^>o (in 
the figure (f> > 4>q, metal is deposited), the barriers are 
modified. To first order in (f> — 4>q, 

AG,. ed (0) = AG red (^o) - a red e((j) - (f> ), 
AG O2; (0) = AG ox (<t>o)+a ox e(c/)-<t) ), (A3) 



FIG. 17: Free enthalpy along the reaction path for the elec- 
tron transfer. Between the two states M° and M + + e, 
there is a transition state J of higher enthalpy. There- 
fore, the barriers AG re( j and AGox determine the reaction 
rate. It exists an equilibrium potential (f>o for which M° and 
M + + e have the same enthalpy, the two barriers are the 
same, AG re d(0o) = AG»i(^o), and there is a balance be- 
tween oxidation and reduction. When a potential cj> 7^ <j>o 
is applied, the relative positions of the states vary, and the 
barriers are modified. In the figure, the solid line is the free 
enthalpy at <f> 7^ <f>o, whereas the two dotted lines are two 
copies of the equilibrium enthalpy profile, shifted to match 
the two local minima in the non-equilibrium profile. The def- 
initions of the modified barriers can be read off; since we have 
a red e((j) — (f> ) + a ox e(<f> - (f>o) = e(4> - (po), Eq. HA41 follows. 
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APPENDIX A: THE BUTLER- VOLMER 
EQUATION 

On the metal-electrolyte interface, the oxido-reduction 
reaction, 



M+(k) + e~(k + a) f± M°(k). (Al) 



with (see Fig. El for an explanation), 



(A4) 



The Butler- Volmer relation then gives the electron trans- 
fer current 



J = Jo exp 



ea red ( 



.) 



kT 



exp 



ea x(4> - <M 
kT 



(A5) 

In our present model, we have supposed in Eq. (|A1(I that 
the reduction of a cation in k is due to a charge transfer 
from the metal site k + a. The link with Eq. (|23|) for the 
reaction rate can then be established. The electrochemi- 
cal potentials are 



A, 



where we have introduced the potential difference 



4> 



k+a-»k 



= ^k+a " Vk " 



eV{E F )' 



(A6) 



(A7) 



The equilibrium (absence of reaction) is obtained when 
$k+a^k = * k +a^k for which 

K + K+* = ^t+E F - e$ k+a ^ k = Mk . (A8) 

Furthermore, a reaction rate of the form 1231) can also be 
written as 



is characterized by two rates k ox and k red . In our sim- 
plified model, the reduction of a cation M + located on 
a site k that is nearest neighbor of a surface site, is car- 
ried out by the transfer of an electron coming from a site 
k + a of the electrode: a reduced metallic atom is then 
created at site k on the interface. 

The Butler- Volmer model (Fig. I17|) supposes that 
there exists a potential difference 4> between the elec- 
trolyte, where the cations M + can be reduced, and the 
metal of the electrode, where a metallic atom M° can be 
oxidized. In between, there exists an activation barrier 



c k , k + a 



with 



k,k+a 



exp 



exp a red - 



y k,k+a — w k,k+a ex P ' 



kT 



Mk + /4+a - A 

kT 



(^k + Mk4 



(A9) 



£*redMk 



kT 



(A10) 
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With the help of Eq. i|A6|l , we obtain 

eaw^k+a^k - $ k +a^k) 



Ck.kH 



y k,^ 



exp 



kT 



exp 



e« OI (*k+a^k-^ +a ^ k ) 



kT 



which has the form of Eq. QA5JI . In the above expressions, 
for a square or a simple cubic lattice, with intersite dis- 
tance a, the transfer current density is 



a d ^k.k+aexp 



kT 



y k,kH 



(A12) 



where the last expression is valid close to equilibrium. In 
that case, and with the help of Eq. (Pljl for w k k+a , the 
(AH) constant jo of the Butler- Volmer law can be identified as 
jo = a d ~ 1 w*e tl l kT , where /jP is the equilibrium chemical 
potential for the metal species. Note that with Eq. I|24|l . 
our expression is only an approximation to the Butlcr- 
Volmer law, valid close to equilibrium; to get a complete 
correspondence with the Butler- Volmer model, a depen- 
dence of w* on the electrochemical potentials needs to be 
introduced. 
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